function wmvmr = vmr_MarkWhiteMatter(hfile, opts)
% VMR::MarkWhiteMatter  - mark white matter (e.g. for IC)
%
% FORMAT:       wmvmr = vmr.MarkWhiteMatter([opts])
%
% Input fields:
%
%       opts        optional settings
%        .gsmoothc  gradient smoothing cycles (default: 6)
%        .smoothc   smoothing cycles (default: 5)
%        .smoothr   smoothing radius (default: 5)
%

% Version:  v0.7c
% Build:    7120608
% Date:     Dec-06 2007, 8:20 AM CET
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'gsmoothc')
    opts.gsmoothc = 6;
end
if ~isfield(opts, 'smoothc')
    opts.smoothc = 5;
end
if ~isfield(opts, 'smoothr')
    opts.smoothr = 5;
end

% copy object
vc = aft_CopyObject(hfile);

% get content of copy
bc = bvqxfile_getcont(vc.L);
if ~isempty(bc.VMRData16)
    bc = bc.VMRData16;
else
    bc = bc.VMRData;
end
bm = mean(bc(:));

% perform number of gradient smoothings
for sc = 1:opts.gsmoothc
    bc = gradsmooth(bc);
end
bd = bc * (bm / mean(bc(:)));

% set content to temporary VMR
bc = bvqxfile_getcont(vc.L);
if ~isempty(bc.VMRData16)
    bc.VMRData16 = uint16(bd);
else
    bc.VMRData = uint8(bd);
end
bvqxfile_setcont(vc.L, bc);

% create gradient VMR
g = vmr_GradientVMR(vc);

% get gradient data
gc = bvqxfile_getcont(g.L);
gc = gc.VMRData;

% delete temporary VMRs
bvqxfile_clear(g.L);
bvqxfile_clear(vc.L);

% mark white matter
wm = (bd > mean(bd(:))) & (gc < mean(gc(:)));

% create WM VMR
wmvmr = BVQXfile('new:vmr');
wmc = bvqxfile_getcont(wmvmr.L);
wmc.VMRData = uint8([]);
wmc.VMRData(size(wm, 1), size(wm, 2), size(wm, 3)) = 0;
wmc.VMRData(wm) = 240;
bvqxfile_setcont(wmvmr.L, wmc);
